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Reconstructing complex networks 
from time series 



Zoran Levnajic 
Faculty of Information Studies, Novo mesto, Slovenia 



Abstract. 

Novel method of reconstructing the topology of dy- 
namical networks from time series is proposed. By ex- 
amining the variable-derivative correlation of the net- 
work node pairs, we derive a simple equation yielding 
the network adjacency matrix. Our key assumption is 
that the intra-network interaction functions are known. 
We illustrate the method on a simple example, and dis- 
cuss the dependence of the reconstruction on the dy- 
namical properties of time series. Our method is appli- 
cable to any weighted or directed network, in principle 
allowing for precision to be estimated. 



1 Introduction 

Complex systems are ubiquitous in nature. On all scales 
from genes to societies, various real systems are composed 



work dynamics. Non-invasive reconstruction methods fo- 
cus on investigating the observable network outputs, such 
as time series of the node's dynamics, or similar measur- 
able quantities describing the internal network states [1011 1112113] . 
The relevance of the non-invasive approach is increasingly 
recognized, particularly due its suitability for detecting 
links in biological networks [12113] , However, by not need- 
ing to interfere with the network dynamics, non-invasive 
approaches are typically limited to studying equilibrium 
network behavior, which often contains only a small amount 
of information on the system. Alternatively, reconstruc- 
tion methods also rely on other mathematical techniques, 
such as control theory [14] . compressive sensing [15] . and 
even treat non-equilibrium scenarios [16] . 

In this contribution, we propose a novel non-invasive 
network reconstruction method. Classical studies usually 
involve correlations among the dynamic variables, i.e. net- 
work nodes [T7] . On the other hand, typical models of net- 
work dynamics rely on first-order differential equations (in 
genetic interactions for example, it is known that one pro- 
tein's concentration determines the increase or decrease 
of another protein's concentration [IB])- Inspired by this, 
we examine the correlation between the variable of one 
node, and the derivative of another node. Our central as- 
sumption is the precise knowledge of the functional forms 
of the intra-network interactions. As we show, depend- 
ing on the quantity of network information contained in 
the empirical data, our method can give very precise re- 



of many units which collectively perform complicated tasks poults even for very short time series, thus being suitable 



In the recent years, the framework of complex networks 
became recognized as an excellent formalism for studying 
complex systems. By modeling units as nodes and their 
interactions as links [2], graph analysis methods entered 
physics, biology, engineering and even sociology [3]. This 
allowed for a variety of real and artificial complex systems 
to be extensively examined, typically via computational 
modeling [3]. Crucial aspect of a complex network is its 
structure, i.e. the topology of connections among its nodes. 
Properties of network structure dictate its global behavior, 
and are key to understanding the network's functioning 
and potentials for its control. For phase-repulsive oscil- 
lator models, profound intcrtwinement between network 
structure and network dynamics was recently shown [5|. 

Since the structure of many natural networks is only 
partially known, it is of central interest to develop mcth- 



for actual experimental situations. Apart from being non- 
invasive, our method is conceptually very simple and easy 
to numerically implement. Based on similar hypothesis, 
Shandilya and Timmc recently proposed another recon- 
struction method [10] , In contrast to their result, our method 
avoids solving the overdetermined linear system, and in 
principle allows for the reconstruction error to be esti- 
mated. 

The paper is organized as follows: after exposing the 
reconstruction method in next Section, we illustrate its 
implementation via simple example in Section^ We present 
the possible extensions of our method, and discuss its lim- 
itations in Section @] 



2 The Reconstruction Method 



ods for reconstructing the network topology from the avail- 
able empirical information. Various experimental techniques We consider a complex system composed of N interact- 
in this directions are already in use, specially in the con- ing units, which we represent as a network with N nodes, 



text of gene regulation networks [H| • In addition, a range of 
mathematical results is available [7]. Recently, the topol- 
ogy of a social network was inferred using mobile phone 
data [5]. Theoretical reconstruction methods usually rely 
on examining computational network models, and are gen- 
erally divided into two classes. Invasive methods involve 
interfering with the network dynamics via controllable 
perturbation, which allows for structural data to be easily 
extracted [5] . These methods generally give very good re- 
sults, specially when the usage of perturbation is not lim- 
ited in strength and frequency. However, it is often unprac- 
tical or even impossible to interact with the on-going net- 



whosc links model the node interactions. Each node is 
assigned a dynamical state defined by the real variable 
Xi = Xi(t), where i = 1,,..N. We assume to be in the pos- 
session of empirically obtained trajectories Xi(t m ). They 
describe the network's discrete-time dynamical evolution 
over a certain time interval. The available data consists of 
N sequences, each containing L values Xi(t\), . . . x»(ti). 
The measurements of Xi are separated by the observation 
interval St = — tm, which defines the resolution of 

the time series (sampling frequency). Time interval St is 
uniform and assumed smaller than the characteristic dy- 
namical time scale. 
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Wc further assume the time-evolution of the node i to 
be governed by: 



N 



(1) 



where we describe the local dynamics via function /l(x), 
and the network (inter-node) coupling by the function fc- 
The network structure is encoded in the adjacency matrix 
Aij , whose element ij specifies the strength with which the 
node i acts on the node j. The complete dynamics of the 
node i is thus the cumulative effect of the local dynamics 
and the contributions from its networks neighbors that 
come with different strengths. Finally, we also assume that 
both interaction functions Jl and fc are precisely known. 

We seek to reconstruct the network's adjacency ma- 
trix under the named assumptions - by having the 
"fingerprint" of the network behavior (i.e. time series of 
the node trajectories), we attempt to reveal its structure. 
The above assumptions on which we build our theory are 
realistic. Many natural systems are modeled using EqfTJ 
examples include gene regulation and neural interactions, 
for which the interaction functions are widely investigated, 
and do not vary with the network links. Modern experi- 
mental techniques allow for high-resolution measurements 
of quantities such as gene expression, although thus ob- 
tained time series are typically short. 

Inter-depcndcncc between the dynamical quantities is 
usually quantified as their correlation |17j . Inspired by 
this, we investigate the correlation between a node vari- 
able (xi) and the derivative of another node variable (xj). 
We define the following matrices: 



B, 



(XifL(Xj)) , 
(XifciXj)) , 



(2) 



where (•} denotes the time-average of a dynamical quan- 
tity (i.e., average over the recorded time-evolution) (h) = 

j- J^m=i M*m)* This allows for Eq|I] to be re-written in 
the matrix form: 



To discern from the original adjacency matrix A^, we 
term the reconstructed adjacency matrix i?^, and quan- 
tify the matrix reconstruction error as follows: 



A 



M 2 



V A 2 



Natural test to make for each obtained Rij is to quan- 
tify how well does it reproduce the original empirical data 
Xi(t m ). To achieve this, we apply the following procedure: 
start the run from Xi(t\) for all nodes, and run the dynam- 
ics using reconstructed matrix for the time interval 5t , 
i.e. until the time i2- Denote thus obtained values l/ifa), 
re-start the run from Xifa) running until £3, accordingly 
obtaining yi(t 3 ), and so on. The discrepancy that the time 
series yi(t m ) show in comparison to Xi(t m ) is the simplest 
measure of the reconstruction precision. We name it tra- 
jectory error At and define as follows: 



An 



N ^ 



Tmkt*".) -Vi{t m )} 2 



This way we measure point-by-point exactness of the re- 
constructed trajectory, which quantifies how well does it 
conform to the actual data. As we show in what follows, 
two errors are in general related, meaning that small At 
suggests small Aa- 



3 Results 

We test our reconstruction method using a simple illustra- 
tive example. A network with N = 6 nodes is constructed 
by placing L — 17 directed links between randomly cho- 
sen pairs of nodes, while requiring the resulting network to 
be connected. Links are weighted with positive and nega- 
tive weights, uniformly selected at random from [—10, 10]. 
The resulting network, illustrated in FigfTJ is thus both 
weighted and directed. Its adjacency matrix Ay is shown 
in Fig[3^,. The dynamics is defined on the network via 



At 



(Bik — Cik) ■ E 



(3) 



which is our main network reconstruction equation. To 
obtain a more stable function and derivative estimates, 
we introduce a new set of time points |10j : 



1, 



,L-1, 



so that ±i(T m ) = [xi(t m+1 ) - Xi(t m )]/6 t and accordingly 
fL,c{T m ) = [fL,c(t m +i) + fL,c(t m )}/2. We will rely on 
this calculation scheme for the implementation of our the- 
ory through EqI2] In principle, our method is applicable 
to any network, and the reconstruction is precisely cor- 
rect in the limit of very long time series. However, since 
the empirical data are not only finite, but typically very 
short, our method will in general yield an approximate 
reconstruction. 




Fig. 1. Graphical representation of the studied net- 
work. Link thickness illustrates the interaction strength. 
Red/blue link colors (light /dark shades) indicate posi- 
tive/negative inter-node interactions. Nodes are numbered 
in accordance with the adjacency matrix shown in FigJSh-. 
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Hanscl-Sompolinsky model [TH] by putting Jl = ~x and 
fc = tanh x in EqJT] The complete dynamics on network 
reads: 



6 



(4) 



For each node we randomly select an initial condition from 
[— 1, 1], and numerically integrate Eq@]from time t = to 
t = 3. During the run, we store 15 values for each Xi, 
equally spaced in time, starting with Xi(t\ = 0). Thus 
obtained time series for all nodes are shown in FigO 




Fig. 2. Time series for all 6 nodes for network FigfTJ ob- 
tained for the first set of initial conditions. 
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Fig. 3. 

trix Rij reconstructed from time series in Fig|2] in (b). 
Colorbar (shade) indicates the link weights. Link-by-link 
comparison of (circles) and Rij (crosses) in (c). Matrix 
error Aa = 0.18, trajectory error At — 0.038. 




We assume now these time series to be obtained from 
an "external" source, such as an experimental measure- 
ment, and use them to reconstruct the network's adja- 
cency matrix as discussed in the previous Section. To this 
end, we numerically compute the derivatives Xi{r m ) and 
the matrices Bij, C'ij and Eij, in order to obtain the re- 
constructed adjacency matrix Rij via Eq[31 The result is 
shown in Fig|3]- the original Aij in (a), the reconstructed 
Rij in (b), and link- by-link comparison of A^ and Rij 
in (c). The reconstructed matrix Rij approximates the 
original Aij rather well for all values of link weights. The 
matrix error is Aa = 0.18, and the trajectory error is 
At = 0.038, indicating a good reconstruction precision. 

We now run another simulation of our dynamical sys- 
tem EqUwith the same underlying network, but this time 
starting from a different set of initial conditions. The new 
time series of equal size and resolution is obtained and 
shown in FigHJ We use them to reconstruct our network 
again, and compare two reconstructions of the same net- 
work, obtained via two sets of time series. The new results 
are shown in FigJSJ and organized in analogy with the 
Fig[3J The new Rij has the matrix error Aa = 0.56 and 
the trajectory error of At = 0.050, which is considerably 
worse than in the previous example. The reconstruction 
errors are bigger on most of the links, as it can be clearly 
seen by comparing Fig(3J: and FigEt- 

Despite that both time series are originating from the 
same dynamical network, two reconstructions are differ- 
ent, one of which is closer to the actual network. This 



Fig. 4. Time series for all 6 nodes for network FigJTJ ob- 
tained for the second set of initial conditions. 



shows that besides the length and resolution of the time 
series, the reconstruction precision crucially depends on 
the "quality" of the time series as well, i.e. on the quan- 
tity of network information contained in them. Time series 
showing generic patterns (e.g. periodic oscillations) are 
more easy to reproduce, than the time series displaying 
more peculiar dynamics (e.g. strong chaos). The former 
dynamical data contains less extractable network informa- 
tion than the latter data. In the above two cases, the first 
time series contain more network information than the 
second, although this cannot be directly seen from Figj2] 
and Fig|2J which explains the difference in precision. 

On the other hand, it is generally difficult to direct 
the experimental measurements towards maximizing the 
available network information. Is there a way to optimize 
the reconstruction precision for any given time series? The 
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(c) 

Fig. 5. Original adjacency matrix Ay in (a). The ma- 
trix Rij reconstructed from time series in Fig f?] in (b). 
Colorbar (shade) indicates the link weights. Link-by-link 
comparison of Ay (circles) and i?y (crosses) in (c). Matrix 
error Aa = 0.56, trajectory error At = 0.050. 



final Section of this paper is devoted to answering this 
question. 



4 Discussion and Conclusions 



The proposed reconstruction method applies to any net- 
work whose internal interactions are described by EqflJ 
and whose interaction functions and fc are known. 
However, the final reconstruction precision depends on a 
number of factors: (i) length and resolution of time series, 
also related to the precision of derivative estimates; (ii) 
quantity of network information contained in the empiri- 
cal data, which can be seen as the reproducibility of time 
series, or coverage of the dynamical phase space with data; 
(Hi) noise in the systems, which can manifest itself either 
as an error in the measurement of Xi(t m ) or imprecision 
in fi, and fc', (iv) invertibility of the matrix Sy, which 
becomes unstable for det Eij « 0; and finally, (v) prop- 
erties of the network itself - some networks are more re- 
constructable than others, due to exhibiting qualitatively 
different collective dynamics. In a concrete reconstruction 
problem, it is difficult to isolate how much each factor con- 
tributes to A a- Instead of quantifying this, we propose a 
generalization of our method, done towards improving and 
controlling the reconstruction precision, regardless of all 
named factors. 

Our method is based on calculating the correlations 
between the variable Xi and other terms, as defined in 
Eq[21 More generally, we can replace x\ by g(xi), where g 
is an arbitrary function, without changing the main result. 



EqJ5]now becomes: 

4 S) = (9(xi)xj) , 

Cf = {g{xi)h{xj)) , (5) 
E\f = (gixMxj)) , 

where notation B$ indicates that the matrix -By was 
calculated via Eqj5| using a pre-defined function g. EqJ21 
which now reads: 

!® = {B$-C®)-E®-\ (6) 

still holds for any choice of function g. As before, in the 
limit of very long time series, the reconstruction is pre- 
cisely correct for any choice of g. In realistic scenarios 
involving (very) short time series, the reconstruction pre- 
cision, except being finite, will depend on g. Namely, two 
different g-s will in general yield two different -s, each 
with its own precision. This means that g plays the role of 
a tunable parameter, which can be used to find the best 
reconstruction. Considering various choices of g, one can 
compute R[f for each one of them, and define as the best 
that Rrp whose reconstructed dynamics shows minimal 

AjP . This way, we can manipulate the reconstruction for 
any given time series towards finding the optimal g, and 
hence, the best Kf- . Such Rf} will extract maximal us- 
able network information hidden in the time series, and 
improve the simple reconstruction obtained for g{x) = x. 

Moreover, the variations of R^' with g are related to the 
reconstruction precision. For a reliable reconstruction, the 
obtained R\y will not strongly depend on changes of g. A 
bad reconstruction will be recognized by a drastic depen- 
dence of i?y on g. Note also, that the functional prop- 
erties of g itself are irrelevant - the only role of g is the 
computation of Rfj ' . 

To illustrate the implementation of our generalized 
method, we examine again the second time series shown 
in FigQ] ("lower quality" ones). For simplicity, we choose 
the set of functions g(x) = x n , where for n we take inte- 
gers between -20 and 20 (except 0). The network is recon- 
structed using Eq[5] for each such g, and the correspond- 
ing A^a and A^P are calculated. The results are shown in 
Figl5J where each R[f is represented as a point through 

its A^ and as coordinates. Clearly, different choices 
of g lead to very different reconstruction precisions. There 
is a visible correlation between the two errors, suggest- 
ing that smaller Aj) , on average, leads to a smaller Ajp . 
The best reconstruction, in terms of matrix error A^\ 
is found for g(x) = x~ 6 [A^ = 0.033, A^ ] = 0.10). 
However, the reconstruction displaying minimal Aif' is 
obtained for g(x) = x 19 (A^ ] = 0.020, A { f = 0.11). De- 
spite the minimal A^ not coinciding with the minimal 

A^\ their values are still relatively close. Note that, as 
indicated in FigEl both of these results are much better 
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?(x)=x 




g(x)=x~ 



(g) 



Fig. 6. Reconstructions using many functions g. Each 
point is specified by A)p and A^' corresponding to i2y 
obtained via one of 17(2;) = x n for n integer between -20 
and 20 except 0. Points obtained for g(x) = x, g(x) = x 19 



and g(x) 



are indicated. 



than what initially obtained for g{x) — x. In addition, 
these results are also better than the one found for the 
first time series from Fig J5] ("better quality" ones). This 
indicates, that despite lower network information content 
in the time series from Fig|4j we can considerably im- 
prove the reconstruction precision by adequately tuning 
the function g. We show the reconstruction for g(x) = x 
in FiglXl in analogy with Fig[3]and Fig|5j 
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Fig. 7. Original adjacency matrix in (a). The matrix 
Rij reconstructed from time series in Fig [?] using the func- 
tion g(x) = x 19 in (b). Colorbar (shade) indicates the link 
weights. Link-by-link comparison of A^j (circles) and Rij 

(crosses) in (c). Matrix error A^' = 0.11, trajectory error 

A^ } = 0.020. 



The role of g is to compensate for errors in the recon- 
struction, occurring for the reasons discussed at the begin- 
ning of this Section. Of course, by considering more (lin- 
early independent) g-s, one can improve the precision in 

both A^ and Aj). However, the question of selecting the 
optimal g which extracts all the network information con- 
tained in the time series remains open. Most straightfor- 
wardly, one can search for such g via Monte Carlo method, 
using randomly chosen functions, or through techniques 
such as evolutionary optimization algorithms |20| . On the 
other hand, in the actual reconstruction problem we can 
only measure Ajl\ From Fig(5J this leads to g(x) = x 19 
as our "best guess" for R[f- Here, one can seek to esti- 



mate A*jp by examining the changes of such R^' for small 
variations around g(x) = x 19 . 

A significant limitation of our method is the ubiquity 
of noise. Simple estimation of derivatives is very sensi- 
tive to noise in x%. However, by using the appropriate 
data smoothing techniques [3T], one might still construct 
a reasonable approximation. Noise can also be present in 
form of stochastic time-evolution of our system (additive 
noise term in Eq[T]). Here, if the noise intensity is known, 
one could modify Eq(5]to include this term as well. Our 
strongest hypothesis is the precise knowledge of the in- 
teraction functions Jl and fc- Despite the availability of 
good mathematical models for many natural interactions, 
lifting this assumption would greatly enhance the general- 
ity of our theory. When approximate functional forms are 
known, interaction functions can be expanded in series, fa- 
cilitating their reconstruction. This would mean that for 
each g, we obtain not just it|? , but also and 
This leads to a possibility of obtaining many different net- 
works, all reproducing empirical data equally well, but in 
combination with different interaction functions. 

In conclusion, we presented a novel method of recon- 
structing the topology of a general dynamical network 
from the empirical time series. By investigating the variable- 
derivative correlations, our method reduces to a simple 
matrix equation, as often studied in linear systems the- 
ory [22] . Notably, our method works reasonably well even 
for very short time series, compared to the system size (we 
reconstructed 6x6 matrix from 6 x 15 data points). The 
problem of network reconstruction is similar to the prob- 
lem of designing a network with a prescribed dynamics. 
One can in principle use our method to design a network 
that displays given time series, by specifying the error 

Aj), which here plays the role of tolerance. Of course, the 
key difference between network reconstruction and net- 
work design, is the interpretation of the solutions. In case 
of many different networks that satisfying A^ = 0, any 
of them is the solution of the design problem. In recon- 
struction theory however, one faces the issue of determin- 
ing which one of them is behind the observed dynamics. 
We also note that our key assumption is the mathemat- 
ical form of network interactions given by EqlT] While a 
similar theory could be developed for any known form of 
Eq|TJ the problem arises for networks whose interaction 
form is not (precisely) known, which in often the case in 



?(<?) 
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real (e.g. biological) networks. Finally, the invasive ver- 
sion of our theory would involve time series measured im- 
mediately after an external perturbation. Such additional 
assumption would greatly improve our method, since tran- 
sient dynamics, similarly to chaotic dynamics, contains far 
more extractable network information than steady-state 
dynamics. 
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